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ABSTRACT 



The rationale and numerical technique of embedding an oceanic 
bulk mixed layer model with a multi-level primitive equation model 
is presented. In addition to the usual prognostic variables that 
exist in a multi-level primitive equation model, the embedded 
model predicts the depth of the well mixed layer as well as the 
jumps in temperature and velocity that occur at the base of that 
layer. The depth of the mixed layer need not coincide with any of 
the fixed model levels used in the primitive equations calcula- 
tions. 

In addition to advective changes, the mixed layer can deepen 
by entrainment and it can reform at a shallower depth in the ab- 
sence of entrainment. When the mixed layer reforms at a shallower 
depth, the vertical profile of temperature below the new, shal- 
lower mixed layer is adjusted to fit the fixed- level structure 
used in the primitive equations calculations using a method vhich 
conserves heat, momentum and potential energy. Finally, a dynamic 
stability condition, which includes a consideration of both the 
vertical current shear and the vertical temperature gradient, is 
introduced in place of the traditional "convective adjustment". 

A two-dimensional version of the model is used to test the em- 
bedded model formulations and to study the response of the ocean 
to a stationary axisymmetric hurricane. The model results indi- 
cate a strong interdependence between vertical turbulent mixing 
and advection of heat. 
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1. INTRODUCTION 



Ttie value of numerical models for the study of ocean dynamics 
has become widely recognized in recent years (National Academy of 
Sciences, 1975; Kraus, 1977; WMO/ICSU, 1977) . However, no single 
model has been developed that can resolve the full range of space 
and time scales commonly observed in the ocean. Thus, the charac- 
teristics of each nimerical model differ considerably and they are 
usually tailored to the space and time scales of the oceanic phe- 
nomena being studied. The time scales to which various models have 
been applied range from decades for the deep ocean circulation to 
seasonal, synoptic and diurnal for near-surface phenomena. 

Numerical models of the ocean can be separated into two broad 
categories. These are the basin-wide or global ocean models v\hich 
attempt to simulate the ocean circulation on seasonal or longer 
time scales, and the one-dimensional models which attempt to simu- 
late the locally forced response at seasonal or shorter time 
scales. The former type of model is referred to as an ocean gen- 
eral circulation model (CGCM) and the latter type, vhich usually 
focuses on vertical mixing processes, is referred to as a mixed 
layer model (MLM) . The most recent comprehensive reviews of OGCM's 
and MLM's are those by Holland (1979), Robinson, Harrison and 
Haidvogel (1979), Bryan (1979), Haney (1979) and Garwood (1979). 

The present study, vhich addresses the problem of embedding a 
one-dimensional mixed layer model into a numerical model of the 
general circulation, was motivated by several factors. First, the 
solution of this problem was recently set as an important goal by 
an international group of scientists studying ocean models and 
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t±ieir relation to climate (WMO/ICSU, 1977). In addition, there are 
a wide variety of important ocean problems involving intermediate 
and mesoscale motions vAiich could be addressed with a realistical- 
ly embedded mixed- layer ocean circulation model. These problems 
include studies of open ocean fronts, coastal and equatorial up- 
welling and the ocean's response to intense atmospheric storms. It 
is also anticipated that an embedded mixed layer ocean-circulation 
model will eventually have a number of important applications in 
the problem of operational ocean prediction. 

In this paper, the rationale and numerical technique of embed- 
ding a model of the mixed layer into an ocean circulation model is 
presented. Hie resulting embedded model has been tested in a 
variety of ocean studies. The simulation of the ocean's response 
to hurricane forcing is presented here because the large horizon- 
tal variations in advection and mixing in such a case (Madala and 
Piacsek, 1977) provide a severe test of the embedded model. Model 
simulations of the oceanic response to less extreme forcing, which 
illustrate more subtle aspects of the embedding, will be published 
later. Our purpose is not to intercompare possible mixing para- 
meterizations, but rather to describe the embedding of a particu- 
lar bulk mixed layer model within a multi-level ocean circulation 
model . 

2. MCOEL EQUATIONS AND BOUNDARY CONDITIONS 

In this section, the equations and boundary conditions for the 
embedded mixed-layer ocean circulation model, as applied to the 
hurricane forcing problem, are presented. Details of the numerical 
formulation, however, are given in the Appendix. 
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The ocean is taken to be hydrostatic and incompressible and 
the density is a linear function of temperature alone. For compu- 
tational economy and because hurricane forcing is nearly axisymme- 
tric, the ocean's response is assumed to be axisymmetric. In addi- 
tion, the Coriolis parameter is assuned to be constant. Ibe 
governing equations are written in radial coordinates as 
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v^ere u, v and w are the radial, tangential and vertical compo- 
nents of velocity, respectively; p is pressure; T is temperature; 
p is density; Aj^ and are coefficients of horizontal eddy vis- 
cosity and conductivity, respectively; a is the coefficient of 
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thennal expansion; and is the density of water at the reference 
temperature, T^. The values of constants used in the hurricane 
forcing experiments are given in Table 1. The terms (w' u' ) , 

(w' v' ) and (w'T* ) represent ensemble average turbulent fluxes and 
are described telow. In (2.1)- (2.3), the advection terms have 
been written in flux form using the continuity equation (2.4). 

The above model equations are applied to an oceanic region ex- 
tending from the surface to a depth of H = 964 m and from the 
storm center (r = 0) to a radial distance of R = 450 km. There 
are no fluxes of mass, momentum or heat across the boundaries ex- 
cept at the surface v^ere a steady wind stress and upward heat 
flux are prescribed. The specific conditions are 




(2.7) 
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The forcing functions Q, t and x are shown in Fig. 1. The 

r 9 

stress is calculated frcxn a prescribed wind, the magnitude of 

viiich increases linearly with r from the eye wall of the storm 

(r = 4.5 km) to the radius of maximum winds (r = 45 km); de- 
- 1/2 

creases like r from there to 360 km; and then decreases 

linearly to zero at the boundary (450 km) . This wind profile 

results in the wind stress components shown in Fig . 1 and a 

curl, 3(t r)/rar, vhich is zero inside the eye of the storm; 

0 

positive and increasing linearly from r = 4.5 km to 45 km; zero 
from 45 km to 360 km and negative from r = 360 km to the bound- 
ary. The maximun tangential stress of nearly 36 dPa corresponds to 
a maximum surface wind speed of 50 m s The small inward radial 
stress component, x^. < 0, results from a cross isobaric surface 
inflow angle of 20 degrees. The effect of other choices of wind 
stress profiles will be reported elsev^i^ere. 

Putting w = 0 at the top of the ocean, as in (2.8) , is known 
as the rigid-lid approximation and it requires the divergence of 
the vertically averaged motion to be zero. In the present model, 
we simplify the dynamics even further by assuming that the ver- 
tically averaged motion itself is zero. This simplification is 
justified in the present case because the shear currents spun up 
by the surface stress are much greater than the dynamically in- 
duced barotropic currents. Subtracting the vertical average of 
(2.1) and (2.2) from (2.1) and (2.2) respectively gives the pre- 
diction equations for the vertical shear currents 
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Viiiere ( ) denotes a vertical average over the total ocean depth 

H, p is p-p, and (u,v) are now the shear currents having zero ver- 
tical average. By integrating (2.5) and using (2.6), p is obtained 
as a function of T alone, thus 
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The turbulence closure method described below requires a prog- 
nostic equation for the depth of the well-mixed layer, denoted by 
h. Ibis equation is derived by integrating the continuity equation 
(2.4) over the mixed layer, and applying the boundary condition 
w = 0 at z = 0. Thus, 
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The right hand side of (2.15), denoted w^, is an entrainment ve- 
locity \^hich is related to the net acquisition of mass in the 
mixed layer. Its parameterization is described below. The six 
equations (2.11), (2.12), (2.3), (2.13), (2.4) and (2.15), along 
with the boundary conditions (2.7)- (2. 10) represent a closed sys- 
tem of equations for the six dependent variables u,v,T,p,w and h 
provided the turbulence quantities, (w'u')/ etc., are para- 
meterized. The values of the turbulent fluxes of heat and momentum 
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at the top and bottom of the model ocean are given by the boundary 
conditions (2.9) and (2.10). The parameterization of the turbulent 
fluxes wi thin the water colunn is given in the following sections. 
3. VERTICAL DIFFUSION 



The classical solution to the problem of parameterizing the 
vertical turbulent fluxes of heat and momentum is the posing of 
eddy viscosity and eddy conductivity formulations such that 



w*u' = 



w'v' 



Sill 

-Sill 



w'T' = - K„ — 
" ^ 3z 



(3.1) 



This is reasonable if gradient diffusion is physically plausible, 
and if the values of and Kp are known to be approximately con- 
stant. This approximation has been shown to produce realistic 
vertical profiles of temperature in the upper ocean below the 
depth of immediate influence of the stirring action of the atmo- 
spheric forcing. Tlierefore, this parameterization is employed 
at all depths below the ocean surface mixed layer, with and 
both equal to 0.5 an s~^. It is possible that the dynamic insta- 
bility parameterization described later could alone provide the 
desired mixing below the turbulent boundary layer. 

The above parameterization of the vertical fluxes is not sa- 
tisfactory, however, between the ocean surface and the top of the 
stable thermocline. In this region of intense turbulent mixing, 
the appropriate values of and would be two to three orders 
of magnitude larger and time dependent. Even if realistic values 
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of and can be computed from scaling argiraents or from more 
sophisticated turbulence closure assumptions, the most important 
problem is to determine the depth below which the intensity of the 
turbulence becomes insignificant. Ttie position of this disconti- 
nuity is at z = -h, the depth of the bottom of the ocean surface 
turbulent boundary layer or mixed layer. Thus an entrainment 
model is needed to compute h(r,t). 

4. ENTRAINMENT AND MIXED LAYER MODEL FORMULATION 

In the GarwDod (1977) model, the buoyancy flux at the base of 
the turbulent boundary layer is modeled according to 

ag (w'T') , = - < E >/h . (4,1) 

-n 

This expression is derived from the local (z = -h) turbulent 
kinetic energy budget. In the entrainment zone, the buoyant damp- 
ing attributable to (w'T')_^ is assumed to be balanced by the con- 
vergence of turbulent kinetic energy flux. This convergence term 
is modeled as a function of the turbulent kinetic energy compo- 
nents and the vertical distance over which this energy must be 
transported. Here <w'w'>, and <E> are the bulk (mixed layer aver- 
age) values of the vertical component and total turbulent kinetic 
energies respectively. The prognostic equations for these mixed 
layer turbulence variables are derived using the bulk second order 
closure methods of Garwood (1977) and are 

1^ (h<E>/2) = - agh (w'T ' )_^/ (2Ri*) 

- (<E>^^^ + fh) <E> (4.2) 
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1^ (h<w'w*>/2) = agh ((w'T’)_^ - (w*T')^)/2 

+ (<E> -3<v?^>) + fh) <E>/3 , (4.3) 

1/2 

In (4.2) and (4.3), u* = (t/p^) with x the magnitude of the 

* 

surface stress and p the (constant) air density, and Ri = oghAT/ 

2 2 

(au + Av ) is the bulk Richardson number. The quantities AT, Au 
and Av are the temperature and velocity jumps, or differences, 
between the values in the mixed layer (assumed to be vertically 
uniform) and the underlying stable water, itiese bulk turbulent 
kinetic energy equations are solved algebraically by assuming a 
quasi-steady state for the turbulent kinetic energy budgets: 

(h <E>) = 0 ; (h <wV>) = 0 . 



This assumption causes the computed values of <w'w’> and <E> to be 

in piiase with the surface boundary conditions — the wind stress, 

2 

p^u* , and the effective surface buoyancy flux, ag(w'T')^. In a 

study of nonstationarity including a diurnal cycle (Garwood and 
Yun, 1979) , this approximation was found to be sufficiently accu- 
rate for time scales greater than a few hours. Ttie motivation for 
this quasi-steady state assunption is a substantial gain in numer- 
ical efficiency. The time step is limited more by the requirement 
of resolving the diurnal forcing than for accuracy in solving the 
turbulent kinetic energy equation. This assumption also makes the 
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boundary layer model more compatible with the nunerical require- 
ments of the ocean circulation model. 

For some applications, the bulk Richardson number term in 

★ 

(4.2) may also be neglected if Ri » 1. However, this term 
should be retained in situations in Vvhich strong shears might be 
expected across the entrainment zone. Niiler (1975) and DeSzoeke 
and Rhines (1976) demonstrate that this term may be significant at 
the outset of a deepening event, provided the initial mixed layer 
depth is small. 

(>ice the entrainment buoyancy flux is knowi from (4.1), the 
downward fluxes of momentun associated with entrainment at the 
base of the mixed layer, (w'u')_j^ and (w'v')_j^, are given by 

- (v?rr)_^ = w^Au = - (w'T')_j^ ^ (4.1a) 

- (w'v')_j^ = w^Av = - (w'T')_^ , (4.1b) 



v^ere the entrainment velocity, 



W = rfs 



(4.1c) 



The above mixed layer formulation was selected for coupling 
with an ocean circulation model because of its general applica- 
bility and because it has been extensively tested and calibrated 
at several weather ship stations (Garwood and Adamec, 1980) . The 
general applicability of the above model is attributed to two 
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important features: (1) An annual-period cyclical steady state is 
possible because of a dissipative limitation to deepening that is 
dependent upon planetary rotation; (2) Entrainment flux is a func- 
tion of implicitly computed values of both horizontal and vertical 
components of turbulent kinetic energy. As a result, the fraction 
of turbulent kinetic energy that is available for mixing at the 
base of the mixed layer is not a constant but is dependent upon 
both the surface buoyancy flux and the surface friction velocity. 
These desirable properties are obtained using second order closure 
techniques without sacrificing the bulk model's nonerical effi- 
ciency vAiich is necessary for practical application of the final 
coupled system with large two and three dimensional arrays. 
Boundary layer deepening by entrainment 

The entrainment heat flux, (w'T' )_^, is determined algebra- 
ically from (4.1), and the steady state forms of (4.2) and (4.3) 

* 

using given values of u*, (w'T')^, h and Ri . This heat flux is 
then imposed upon the given temperature profile. The new h and 
the new momentum and temperature profiles are unique solutions 
provided; (i) the mixed layer (0 > z > -h) is homogeneous; (ii) 
the value of h is made just large enough to prevent unstable den- 
sity profiles (8T/3Z > 0); and (iii) heat is conserved. Note that 
the depth of the bottom of the mixed layer is not forced to coin- 
cide with any of the prescribed model levels or interfaces between 
levels. As a result, there is added vertical resolution in the 
vicinity of the mixed layer interface. This extra resolution is 
very desirable because the vertical structure of the thermocline 
determines the potential energy of the upper ocean v\hich is 
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important for both the dynamics of the boundary layer and the 
larger scale baroclinic circulation. 

Boundary layer shallowing 

Equally important in a mixed layer model of general applica- 
bility is the prediction of boundary layer shallowing. This phe- 
nomenon occurs vshen there is inadequate vertical turbulent energy 
to transport heat down to the base of the existing mixed layer. 

The (new) depth at v^hich the vertical flux (w'T') vanishes es- 
tablishes the bottom of the turbulent boundary layer and vAiat will 
be the base of the new mixed layer. In the present model, this 
occurs v^hen <w'w'> approaches zero. Then (4.1) no longer applies 
and the steady state forms of (4.2) and (4.3) reduce to the 
following diagnostic equations for the new mixed layer depth 

0 = + fh) <E> (4.2a) 

0 = -agh(^?T^^/2 + - (<S>^^^ + fh) <E>/3 . (4.3a) 

Without the term involving the planetary rotation, or if the down- 
ward surface buoyancy flux is dominant, h is proportional to the 
Obukhov length scale, L = u* /(ag (w'T' ) ^) . On the other hand, if 
(w'T')^ = 0, h is proportional to u^/f. In general, however, the 
depth of the shallowing mixed layer is a function of the two non- 
dimensional parameters, h/L and hf/u*. 

In a numerical model, the formation of a new, shallover mixed 
layer is potentially a problem because the interfacial structure 
of previously created mixed layers can not be exactly preserved 
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with a relatively coarse vertical grid. In the present model, 
this potential problem is avoided using a numerical procedure (see 
Appendix) that conserves heat, momentum and potential energy, 
vviiile minimizing changes in the vertical density structure. 
Although potential energy is conserved during the shallowing pro- 
cess, no attempt is made to exactly conserve mean kinetic energy. 
5. DYNAMIC SIABILITY CONDITION 

In the derivation of the above bulk mixed layer model, it is 
assumed that the turbulent boundary layer is dynamically unstable 
and that the underlying water colunn is dynamically stable as 
measured by the gradient Richardson number, Ri. In the entrain- 
ment zone (-h _< z _< -h-d) , dynamic stability is assumed to be near 
neutral; that is Ri = Ri™ = 1/4. This condition is insufficient 
in itself to determine the entrainment fluxes or the rate of mixed 

layer deepening; however, it does lead to a measure of the thick- 

2 2 

ness of the zone; d = ■*" 4v ) ogAT. Deepening of the mixed 

CH 

layer depends upon the added intermittent vertical flux of 
turbulent kinetic energy from above. 

As mentioned earlier, dynamic instability can sometimes occur 
below the mixed layer, and this will enhance the vertical trans- 
port of both heat and momentum. Therefore, in the present model, 
dynamic stability is enforced throughout the water column by im- 
posing those vertical fluxes of heat and momentum between model 
levels v\hich are necessary to prevent the gradient Richardson num- 
ber from going below the critical value, i.e. Ri(z) > Ri™* This 

— CK 

"dynamical adjustment" is a simple generalization of the "convec- 
tive adjustment" parameterization which is conmonly used in most 
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ocean general circulation models. Details are given in the 
Append ix . 

It should be pointed out here that the above dynamical sta- 
bility condition is different from that of Pollard, et al . (1973) 

and Thompson (1976) . In that model, deepening is dependent upon 

★ 

stability as determined by the bulk Richardson number, Ri . This 
is equivalent to the present requirement that there be a constant 
gradient Richardson number in the entrainment zone only if the 
thickness of the entrainment zone divided by the mixed layer 
depth, i.e. dA>f is a constant. There is no physical justifica- 
tion for such an assumption. 

6. METHOD OF EMBEDDING 

The system of equations for the embedded model includes a 
prognostic determination of the mixed layer depth. The mixed 
layer is defined here as being the upper layer of the ocean where 
turbulent mixing is dominant. The temperature and horizontal 
velocities are required to be independent of depth in this region. 

One of the more important parameters in predicting the evolu- 
tion of the mixed layer is the dynamic stability at its base. 

Since density is assumed to be a function of tanperature alone, 
the temperature junp at the base of the mixed layer is critical. 
The dynamic portion of the model is a level model that predicts 
the average of a quantity between two specified depths. The base 
of the mixed layer, however, will rarely lie exactly on a model 
level. Requiring a sufficiently fine vertical resolution to pre- 
vent errors in the thermal and potential energy balances due to 
the truncation of the mixed layer depth to a model level is 
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unreasonable bot±» in terms of computer storage and computational 
time. The problem addressed here is to provide for a meaningful 
communication between the fixed-level, dynamical part of the 
model, and the mixed layer part of the model. 

The coupling of the dynamical and mixing processes in the 
model may be considered in two stages. First, the changes in the 
upper ocean due to advective and diffusive processes are calcu- 
lated in the dynamical part of the model and put into appropriate 
form for the mixed layer model. Second, the changes due to sur- 
face fluxes and entrainment mixing are calculated in the mixing 
part of the model and trananitted to the dynamical part of the 
model. In both of these stages, a special treatment is required 
for the model level which contains the base of the mixed layer. 

At the start of a given timestep, changes due to advective and 
diffusive processes (including the dynamic adjustment described in 
Section 5) are calculated for all the GCM layers in the dynamical 
part of the model and the prognostic variables (u,v,T and h) are 
stepped forward in time (see the Appendix for details of the time 
differencing scheme) . The special treatment for the layer con- 
taining the new mixed layer depth is now described with the aid of 
Fig. 2 in v^ich the dashed line is the vertical structure used by 
the dynamical part of the model and the solid line is the vertical 
structure used by the mixing part of the model. Both the dynami- 
cal and the mixing part of the model use the same structure at all 
levels v^ich lie entirely within the mixed layer. The variable ^ 
represents u, v, or T. If the base of the mixed layer lies within 
the k th model level, the value of ^ in the mixed layer portion is 
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denoted by the value just below the mixed layer is denoted by 
the integrated average over the model level is denoted by 
while the mixed layer depth is denoted by h. The change due 
to advective and diffusive processes is assigned the value calcu- 
lated for the GCM layer iitmed lately above (i.e. the layer bounded 
by ^k-3/2 base of the mixed layer is 

constrained to always lie at or below the base of the first model 
level. The associated change 6^2 is then calculated from the re- 
quirement that the weighted average of and be equal to the 
change 6? predicted for the layer by the dynamical part of the 
model. Ttius, 



6? = 6?^ + Z 2 6^2 f (5.1) 

v^ere Z^^ is the portion of the layer above the new mixed layer and 
Z^ + Z2 = In the circulation model, the advective and diffu- 

sive changes are based on the average properties in each layer. 
Efov^ever, there can be a large velocity and temperature jump across 
the interface h. The use of a layer mean advecting velocity and 
layer mean temperature in the dynamical part of the model neces- 
sarily introduces an error in the calculation of the change 6^, 
v\hich, by our formulation, is introduced into the change assigned 
to 6^2* The error in 6E will be smallest when the depth h is 
close to one of the interfaces between GCM layers, and it will be 
largest when h is near the center of a GCM layer. An analogous 
but slightly more complex formulation is used whenever the dynam- 
ically predicted 6h causes the base of the mixed layer to move 
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into an adjacent model layer. Before proceeding to the mixing 
portion of the timestep, the values of C throughout the mixed 
layer (including are replaced by the vertical average of c 
over the new mixed layer depth. In the particular hurricane simu- 
lation shown below, the small amount of turbulent kinetic energy 
that would be expended (or released) in performing this vertical 
mixing has been neglected in the turbulent kinetic energy balance. 
This is clearly justified in the present case since the amount of 
turbulent kinetic energy generated by shear stress production is 
so large. 

The second stage of the coupling is the mixing portion of the 
timestep. During this stage, the case with entrainment mixing oc- 
curring in association with mixed layer deepening must be distin- 
guished from the case in w^ich the layer reforms at a shallower 
depth. Consider first the entrainment mode in w^ich the mixed layer 
model (Section 4) predicts an entrainment heat flux, (w'T' ) 

< 0. The mixed layer deepening, 6 h > 0, during the timestep At is 
calculated fron (4.1c) . In this case, a volune of water (per unit 
area) equal to fih and having property is entrained into the 
mixed layer having volume (per unit area) equal to h and property 

This results in a new mixed layer of depth h + 6 h and property 

^1 ^ '^^ 1 ' 



Entrainment mixing produces no change in c below the new mixed 
layer depth. If 5 ^ > ^ 2 ' '^ich is always true for the temperature 
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in the model, then 6?^ ^ 0 . If the entrainment deepening does not 
cause the base of the mixed layer to move into an adjacent model 
layer, then the magnitude of the property junp across the mixed 
layer, | 5^ - ^ 2 ^ > will be decreased. In the numerical model, a 
separate array of ^2 values must be stored to permit the jump 
values, - ^ 2 ' ^ always specified for the mixing process. 

This value, and the mixed layer depth h are the only additional 
variables v^ich must be retained in coupling the mixing model and 
the dynamical model (see Appendix for details) . 

Consider now the case of mixed layer shallowing. In this case 
the diagnostic relations ( 4 . 2 a) and ( 4 . 3 a) are solved for the new 
mixed layer depth. The heat and momentum fluxed across the sur- 
face during the timestep are distributed over the new, shallovser 
depth to produce a new value of 5^^. A new ^2 value is determined 
fron the prior 5^. The next advective timestep operates with the 
new, shallower mixed layer. Some information is lost regarding 
the deeper structure because provision is made for but one active 
mixed layer. To preserve some of the information regarding the 
density junp at the base of a prior mixed layer, the vertical 
density profile is adjusted to conserve both heat and potential 
energy viienever the mixed layer reforms at a shallower depth (see 
Appendix for details) . 

7 . APPLICATION WITH STATIONARY HURRICANE FORCING 

The purpose of this section is to demonstrate the capability 
of the embedded mixed layer-ocean circulation model to simulate 
realistic upper ocean changes in response to hurricane forcing. 
Such strong forcing as that of a hurricane is used because a large 
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horizontal gradient in the advective and mixing response is expected 
and this will provide an extreme test of the model. No attempt is made 
to compare the results of the present embedded model to the results of 
other models utilizing different mixing parameterizations. Such a com- 
parison is beyond the scope of the present paper. 

The classical model simulation of the ocean's response to hurricane 
forcing is that of O'Brien and Reid (1967) using an axisymmetric, one- 
layer model. The extension of this type of model to two horizontal 
dimensions is given by Chang and Anthes (1978) , vhile the response of a 
multi-level circulation model is described by Madala and piacsek 
(1977) . As a result of these and other studies, certain features of 
the dynamical response of the ocean to a stationary hurricane are well 
known. The currents in the upper layers of the ocean tend to circulate 
in the same sense as the hurricane winds. Although the hurricane wind 
field is directed not only cyclonically about the storm but also 
radially inward, the net radial movement of surface layer water is out- 
ward due to the Coriolis effect. A vertical-radial circulation, con- 
sisting of intense upwelling within the radius of maximum hurricane 
winds and weaker downwelling at greater radial distance is produced 
by mass continuity. As a result of this circulation, the ocean tem- 
peratures near the center of the storm decrease throughout the stable 
part of the water column that lies below the mixed layer. The ocean 
surface temperatures are also affected if the upwelling water reaches 
the surface. 

There is also an important mixing effect due to the hurricane force 
winds, as indicated by the model studies of Elsberry, Fraim and 
Trapnell (1976) and Chang and Anthes (1978) . Mechanical and convective 
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generation of turbulent kinetic energy in the upper ocean tend to 
achieve maximun values at the radius of maximum winds. As a result, 
the relatively warm surface water is cooled and the coler water imme- 
diately below the initial mixed layer depth is warmed as it is en- 
trained into and is mixed with the surface water. The combination of 
the advective effects and mixing effects described above is clearly 
quite complex and requires high vertical and horizontal resolution to 
simulate. 

lb simplify the testing of the present model, the forcing due to a 
stationary, axisynroetric hurricane is used. A 4.5 km horizontal grid 
spacing is utilized in order to have 10 grid points within the radius 
of maximum winds which is located at 45 km (Fig. 1) . The ocean is 
initially at rest with a uniform surface temperature of 30°C and a 
uniform mixed layer depth of 60 m. Ta reduce the initial shock, the 
hurricane forcing is slowly spun up to full intensity over the first 
12 h. 

The evolution predicted at 12 h and 24 h with the embedded mixed 
layer-ocean circulation model is presented in Figs. 3a-d and Figs. 

4a-d, respectively. Because of the explicit mixed layer, the surface 
inputs of heat and momentum are distributed over a deep layer. Thus, 
the maximum v component at 12 h (Fig. 3a) is rather small considering 
the large magnitude of the surface stress near the center (Fig. 1) . At 
12 h there is a deep layer of cyclonically rotating water in the mixed 
layer and a very snail and weak anticyclonic circulation below the 
mixed layer. A region of relatively weak outward flow is predicted 
(Fig. 3b) for the mixed layer in the inner region. The compensating 
inward flow below is a maximum near the top of the thermocline. As is 
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indicated by the packing of the w-isolines at the center of the hurri- 
cane in Fig. 3c, the associated vertical motion at 12 h is concentrated 
in a very narrow zone. The temperature field at 12 h in Fig. 3d is 
beginning to show the effect of the strong upwelling near the center. 
Tbward the center, the isotherms are tilted upward at progressively 
larger angles. The 29°C isotherm has been lifted about 15 m near the 
storm center. 

At 24 h the tangential (v) components within the mixed layer near 
the center are only slightly stronger than at 12 h. A maximum value of 
less than 1 m s~^ is qualitatively in agreement with the largest cur- 
rents measured as hurricane Eloise passed an instrunented oil rig off 
the coast of Louisiana (R>rristall, 1980). Although the strength of 
the hurricane forcing used here is sonevy^at larger than in Eloise, 
there is a compensating effect due to the larger layer depth in the 
model conpared to the pre^loise ocean conditions of between 30 and 45 
m. At larger r, there is a region of weak, anticyclonic flow that is 
produced as the water moves outward vdiile tending to conserve angular 
monentum. The radial extent of the outward flow near the surface in 
Fig. 4b is murh greater than at 12 h. The maximun of this radial com- 
ponent is nearly as large as the corresponding tangential component, so 
that the current direction is approximately 40° relative to a circle, 
or 60° to the right of the surface stress. The compensating inflow oc- 
curs over a deep layer, althotgh the maximum is clearly located in the 
upper layers of the thermocline beneath the radius with maximum stress. 
The primary upwelling at 24 h remains at the center of the hurricane 
and extends to great depths (only the upper layers of the model are 
shown here) . A much weaker descent occurs at large r to compensate for 
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the surface flow, which has reached the outer boundary by 24 h. The 
primary feature in the temperature field at 24 h (Fig. 4d) is the ex- 
pansion of the zone of upwelled water, v\hich has surfaced near the 
center and has been advected outward. 

The advective effect near the center of the hurricane is so domi- 
nant that it is difficult to isolate the mixing effect. Deviations of 
the predicted temperature fields at 12 h and 24 h from the initial 
values are shown in Figs. 5a and 5b, respectively. After 12 h, the 
maximun temperature decrease occurs near the storm center in associa- 
tion with the strong upwelling there, v^hile a region of temperature 
increase occurs near the outer boundary due to downwelling. There is a 
second zone of increased temperatures at the base of the mixed layer, 
v\hich has deepened as a result of both entrainment and vertical advec- 
tion. A similarly shaped region of temperature increase is found some- 
Vihat farther from the center of the 24 h field shown in Fig. 5b. Such 
a region of temperature decrease above the initial mixed layer base and 
increase below it is an indication of entrainment mixing. It is clear 
frcxn the above temperature changes that in the present case of station- 
ary forcing, advective effects dominate the heat budget. If the hurri- 
cane was translating, upwelling at a point would be followed by down- 
welling and the advective effects would have a tendency to cancel. In 
such a case, the mixing effects would be relatively more important in 
producing a permanent change in the thermal structure. 

The explicit treatment of the dynamics of the mixed layer in the 
embedded model results in a temperature field that is similar to obser- 
vations of oceanic response to hurricanes (e.g. Leipper, 1967). Tbe 
purpose of the limited test showi here is only to demonstrate that the 
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embedding procedure proposed here appears to be capable of treating 
both the mixing effects and the advective effects in a miner ically ef- 
ficient and physically realistic manner. No inference should be made 
regarding the viability of different mixing parameter izations. 

8. CONCLUDING REMARKS 

This paper has addressed the problem of developing a dynamical 
ocean circulation model capable of treating both advective and mixing 
effects in response to realistic atmospheric forcing. The solution 
presented here is the embedding of a bulk mixed layer model within a 
multi-level primitive equation model. In addition, the traditional 
"convective adjustment" often used in such primitive equation models 
has been replaced by a more general dynamical stability condition. The 
resulting model is an attempt to combine the advantages of layer models 
and level models into a single model. 

Other approaches to the problem are clearly possible. One alterna- 
tive approach is to use a higher order diffusion model in v^ich the 
vertical eddy viscosity and conductivity are expressed as functions of 
the model variables. This approach has been utilized by Madala and 
piacsek (1977) using the popular formulas of Munk and Anderson (1948) . 
With such an approach, high vertical resolution is necessary in order 
to represent the sharp transition between the highly turbulent surface 
layer and the weakly turbulent, highly stratified region below. 

Another alternative approach is to develop a multi-layer model vhich 
includes turbulent exchange between the layers. A bulk mixed layer 
could be considered as the uppermost layer in such a model . One of the 
results of the present study suggests that, with this approach, there 
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should also be a turbulent transfer of momentum betv^een the layers, an 
effect that has not been generally included in layer models before. 

The example of the oceanic response to a stationary, axisymmetric 
hurricane examined here is an extreme test of the present embedded 
model because both the advective and the mixing effects are very 
strong. Additional studies with the present model are currently in 
progress. These include the ocean's response to a variety of strong 
atmospheric forcing situations and the response of an upper ocean 
density front to local atmospheric forcing. A three dimensional ver- 
sion of the model is being used to study a variety of oceanic meso- 
scale features as well as large-scale interannual variability (anoma- 
lies) in the central North pacific Ocean. These studies suggest that 
the embedded model is capable of simulating the ocean's response to a 
wide range of oceanic and atmospheric forcing conditions. Other tests, 
including comparisons with observational data, are planned to further 
validate the present model. As a result, it now appears that the model 
can serve as a useful tool to increase our understanding of, and 
eventually our ability to predict, a variety of oceanic phenomena that 
involve both advective and mixing processes. 
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APPENDIX 



Space finite differencing 

The numerical scheme makes use of centered space differencing 
on a staggered grid. In the notation v^^ich follows, the letters j 
and k are integers. The vertical index is k and the variables u, 
V, T and p are defined at level k, k=l, ...K, viiile the vertical 
velocity w is defined at the K-1 interfaces between the levels and 
also at the top and bottom of the ocean v^ere it is zero. 

Standard difference notation is utilized in v^ich for example 

1+1/2, j - ■3i+l,j - ‘5i,j "li+l.j ^ 

Let aZ|^, k=l,...K, be the specified thickness of the k th 
level and 



"k^ ’2 



“ - E ^ 1-1 • 



1; 



(A] ) 



be the value of Z at the base of the k level (Z ^^2 ~ 
^K+1/2 ~ ’ Average layer thicknesses are defined by 






(A2) 



and average layer heights are defined by 






1-1. ... K 



(A3) 
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The finite-difference analog of (2.13) is 



Pk = Pk " P ’ 

v\^ere the wavy overbar denotes the finite-difference vertical 
average 



K 



( ) 



4E 






(A5) 
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and 
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Pk 
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k-1 

E -O 



T )1 
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(A6) 



Thus, P|^ is the thermally induced hydrostatic pressure field v>hich 
has a zero vertical average. 

Let Ar be the (constant) horizontal grid size. Then the veloc- 
ity components u and v are defined at the center of the grid in- 
tervals located at the points 



Ar -t- (j-1) Ar ; j-1, ... J . 



(A7) 



The variables T, p and w are defined at the J-1 interfaces and at 
both ends of the domain (r=0 and R) . Tbese points are located at 



r, 0 
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J-P^ 



= R 




(A8) 



30 



In addition, the space differencing scheme makes use of the fol- 
lowing definitions of the grid volumes. 



”>j,k ■ 2 ■'l 
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The finite-difference analog of the continuity equation (2.4) 
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\A^ere 



u . , = u . , r . AZ, 

j,k j k 



“j+>5.k+!-i ‘ “j m.k+b 



(Al]) 



The finite difference analog of the boundary conditions which are 
needed to apply (AlO) , and the equations to follow, on the bound- 
aries are given in the next section. 

The finite-difference analog of the thermal energy equation 
(2.3) and the prognostic equation for the mixed layer depth (2.15) 
are 
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The finite-difference analog of the radial component of the 
equation of motion (2.11) is 
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() L. J , fv J , «C 



" - ■■r^“'>j,k - 'V/"-'>J.k 



4 (I- -I V. ,,/r.)v, - i^(p)j,ktj aVPo 



(A12) 



(A13) 
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(A15) 
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v^ere. 



J * V > J ‘ J * 
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The finite-difference analog of the tangential component of the 
equation of motion (2.12) has a similar form. 

Boundary conditions 

The lateral boundary conditions will only be given for r=0 
since the conditions at r=R are similar. The variables and T,p 
and w are defined at grid points vhich coincide with the lateral 
boundaries vvAiile (u,v) are defined at points located one-half grid 
distance inside the boundary. Thus, the variables '^ 1/2 k' 

^1/2 k '^1/2 k+1/2 located at the origin while r^^, u^^ ^ and 
Vj^ are located one-half grid distance from the origin. 

The zero normal flux condition at r = 0 is maintained by im- 
posing the following antisymmetric condition on the normal flux; 



u 



0,k 



= - u 



l,k ' 



(A17) 



(A16) 
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\nJiile the condition of no lateral eddy flux of heat and momentum 
at r = 0, vAiich is used in (A12) and (A15 ) , is given by 






(ill-), 






(vr)^ - 0 

'^2 > »< 



(A18) 



Equations (A17) and (A18) are the finite-difference analogs of 
(2.7). 

The finite-difference analogs of the boundary conditions on w 
in (2.8) are 

The finite-difference analog of the boundary conditions on the 
surface fluxes in (2.9) and (2.10) are 



‘ "j+>i 



■(u'u'), , » {(1 )./(>) r Ar 

J > '2 r 1 n j 
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and 



2)^^ ' 2 



(w'u'). = (w'v’) , = 0 , 



(A 20) 



(A21) 



respectively, v»iiere the forcing functions Q, t and t are taken 

r (j 

from Fig. 1. 
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Vertical diffusion 



In this section, and the ones to follow, the numerical methods 
used to calculate the turbulent fluxes in (A12) , (A13) and (A15) 
are given. The vertical flux of heat and momentum due to eddy 
diffusion is calculated at the top and bottom (interface) of the 
model level vhich contains the mixed layer depth and at the top 
and bottom of every model level v^ich lies entirely below the 
mixed layer. That part of the turbulent fluxes in (A12) and (A15) 
vshich is due to vertical diffusion is calculated from 






(w'u'O . , ,, - K, . 6 (u) . , , r. Z.r/AZ, ,, 

-^j.k+’i n j,k-ti, j Id-J; 



(A22) 



with a similar equation for (w'v'). 

Boundary layer entrainment 

This section concerns the numerical treatment of entrainment 
by the embedded mixed- layer model. In a time step At, the amount 
of mixed layer deepening and the associated layer temperature 
change depend upon both the prior thermal structure and the en- 
trainment heat flux (w'T' )_^, computed by solving (4.1) simul- 
taneously with (4.2) and (4.3). Following (6.2), heat is verti- 
cally redistributed, and the changes in mixed layer temperature 
and depth are found as illustrated by Fig. A1 in vshich 

i) 5T = (wT^_^ At/h 
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ii) 6h is determined by requiring 
area B = area A = h 5T . 

Only after applying this entrainment mixing is the new mixed layer 
temperature further adjusted to account for the net surface heat 
flux, (w'T')q (although (w'T' helped to determine the value of 
6T above through its appearance in the turbulent kinetic energy 
equation, (4.3)). Once the new mixed layer depth is determined, 
the horizontal components of mixed layer velocity are adjusted so 
as to conserve momentum and preserve vertical homogeneity in the 
mixed layer. 

Boundary layer shallowing 

When the depth of the turbulent boundary layer decreases, it 
leaves behind a remnant mixed layer structure. In a numerical si- 
muilation of this event, unless the base of the older mixed layer 
exactly coincides with an interface between two model levels, sane 
information muist be lost because only one mixed layer depth is re- 
tained. This loss of information is a potentially serious problem 
because such remnant mixed layer structures influence the poten- 
tial energy in the column. In the present model, when the mixed 
layer depth is predicted to decrease, from h to h' say, the tem- 
perature in the model levels v^^^ich lie partly or entirely within 
the mixed layer are adjuisted along with the temperature in the 
model level immediately below the mixed layer depth (Fig. A2) . 

This adjustment is given by 
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(A23) 
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v^ere k is the model level which contained the original mixed 
layer depth before shallowing, T|^ and h are the model temperatures 
and mixed layer depth before shallowing, while Tj^' is the adjusted 
model temperatures after shallowing. It can be shown that the 
above adjustment of temperatures conserves both heat and potential 
energy. 

The procedure is illustrated in Fig. A2. A slight vertical 
redistribution of heat enables the conservation of potential 
energy. In the case when the new turbulent boundary layer, h' , is 
shallower than 



the new mixed layer depth will not be apparent in the temperature 
profile until the turbulent boundary layer is sufficiently warmed 



k-1 




1=1 
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by a net flux of heat and radiation at the surface. In the illus- 
tration, area C + area D = area E. In this purposely exaggerated 
example, the surface temperature is seen to rise a slight amount. 
However, over a cycle of deepening shallowing and deepening again, 
there will be no net change in surface temperature. 

Dynamic stability condition 

This stability condition, which is a simple generalization of 
the familiar "convective adjustment", is imposed by examining the 
value of the Richardson number, Ri, defined at the interface be- 
tween model levels given below. If the inequality in (A26) is not 
satisfied, the temperature and the velocity components in the two 
adjacent levels are "dynamically adjusted" in accordance with the 
following four assumptions: 

(i) heat is conserved, 

(ii) momentum is conserved, 

(iii) dynamic stability is imposed between the two levels by 
requiring that the Richardson number 



Ri- ^ ^ > (Ri)p-, - .25 , 

(Au) + (Av) ' 



(iv) the mixing ratios for heat and momentum are equal, i.e. 
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where the prime denotes the new, adjusted value. In these expressions 



at, Au and Av are the vertical differences of heat and monentum. 



average of the two levels, 0.5(T|^AZj^ + T|^^2^^k+l^^^^k+l/2 
above stability condition is calculated at the velocity points in the 
model, and therefore, because of the staggered grid, the temperatures 
appearing in (A26) and (A27) are twD grid-point averages. Conse- 
quently, a weighted average of the temperature adjustment determined 
from the above conditions is used to compute the adjustments that are 
applied at the tonperature grid points. With some straightforward 
algebraic manipulations it can be shown that the unique solution to 
conditions (i) through (iv) is 



T|^ - while T, u and v represent the weighted vertical 
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If the mixing ratios in (A27) are denoted by m , it can easily be 
shown that 

^ ^ ' 0 < 1 

m = 1 : Ri £ 0 . (A29) 

Complete mixing between two levels, as defined by m=l, occurs 
only if the value of the Richardson number before mixing is less 
than or equal to zero. 

Time integration 

Time integration is performed using the leapfrog scheme with 

an Euler-backward scheme introduced at the start and at every 11 

timesteps thereafter to remove any solution separation in time. 

The forward differencing schane, however, is used for the eddy 

diffusion and turbulent mixing terms. As indicated in Section 6, 

the calculation of the advection and diffusion terms is split from 

the calculation of the turbulent mixing terms. 

Suppose the dependent variables u,v,T and h are known at the 

two consecutive time levels n-1 and n. When using the leapfrog 

schane, the first step is to calculate the pressure p and the ver- 
* 

tical velocity w using (A4) and (AlO) and the values of T, u and 
V at time level n. Next, the dynamical terms on the right hand 
side of (A12), (A13) and (A15) are calculated, also using the 
variables at time level n. These dynamical terms include the hori- 
zontal advection of u,v,T and h; the vertical advection of u,v and 
T; the Coriolis acceleration; the centripetal acceleration; and 
the pressure gradient force. Next, the horizontal and vertical 
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diffusion tenns are calculated using u,v and T at time level n-1 
and are added to the dynamical terms. The coupling procedure 
given by (6.1) is then applied to these accumulated tendencies 
V(t»ich are then multiplied by 2At and added to the value of the 
variable at time level n-1. This partially advances the fields to 
the new time level n+1. The dynamic stability condition (see 
above and in Section 5) is then imposed and the variables within 
the mixed layer are vertically averaged above the new mixed layer 
depth. Finally, the tendencies due to the surface fluxes (A20) 
and turbulent mixing (Section 4) are calculated and the fields are 
advanced to their final value at the time level n+1. 

When the two step Euler-backward integration scheme is used, 
the sequence of calculations is the same as described above in 
both the "forvrard" and the "backward" steps. The only difference 
is that the diffusion terms are no longer lagged in time but are 
calculated from the same fields as the dynamical terms. 
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TABLE 1 



Symbol 


Value 


Description 




2 X 10® cm^ s~^ 


Horizontal eddy viscosity 




7 2-1 

4 X 10 cm s 


Horizontal eddy conductivity 




.5 cm s 


Vertical eddy viscosity 




2 -1 
.5 cm s 


Vertical eddy conductivity 


a 


2 X lO"'* deg"^ 


Thermal expansion coefficient 


T 

o 


5°C 


Reference temperature 


Pq 


-3 

1.0276 gm cm 


Reference density 


f 


6 X lO"® s“^ 


Coriolis parameter 


hr 


4.5 km 


Horizontal grid size 


r 

m 


45 km 


Radius of maximun winds 


R 


450 km 


Domain extent 


H 


964 m 


Depth of ocean 


At 


150 s 


Time step 
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FIGURE CAPTICWS 



Fig. 1. 



Fig. 2. 



Fig. 3 



Fig. 4 



Ftorcing functions for the ocean model, each normalized by 
the value of the function at the radius of maximun winds. 

The tangential and radial stress components have the same 
form (solid line) with maximum values of 35.9 and -12.9 dpa, 
respectively. The surface heat flux (dashed) has a maximun 
value of -840 Wn~^. 

Schematic illustration of the emtiedded mixed layer of depth 
h, located within the k th level in the circulation model. 
See Section 6 for a complete discussion of the embedding 
technique . 

. Predicted oceanic response after 12 h to the hurricane 
forcing in Fig. 1 using the circulation model with the 
embedded mixed layer. Only the upper 150 m is shown, with 
the center of the hurricane (r=0) on the left and each 
larger tick mark on the abscissa corresponding to 45 km. 

(a) Tangential current component (cyclonic flow, dashed) 
with contour interval of 0.1 m s~^ . (b) Radial current 

component (inflow, dashed) with contour interval 0.1 m s“^ . 
(c) Vertical velocity (downvvard flow, dashed) with contour 
interval of 1 m h”^ . (d) Tanperature contours (°C) . 

. Similar to Fig. 3 except at 24 h. 
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Fig. 5. Temperature deviations (contour interval 0.5°C, negative 
values dashed) from the initial state predicted by the 
embedded mixed layer model after (a) 12 h and (b) 24 h. 

Fig. Al. An example of mixed layer deepening and cooling by entrain- 
ment. The original structure prior to entrainment (solid 
line) had a mixed layer depth of 36 m, mixed layer tempera- 
ture of 22.5°C and temperature jump at the base of the layer 
of approximately 2.6°C. The mixed layer deepened by entrain- 
ment to fc>elow 40 m, the mixed layer temperature decreased to 
22.2°C and the temperature jimp at the base of the layer 
decreased to approximately 2.4°C (dashed line) . Heat is 
conserved, so area A = area B. 

Fig. A2. An example of mixed layer shallowing. The original tempera- 
ture structure prior to shallowing (solid line) is the same 
as in Fig. Al. In this example, the mixed layer shallowed 
to a new depth which lies above the k th level, but vhich 
will not be apparent in the temperature structure (dashed 
line) until there is a dov\nward flux of heat at the surface. 
The new temperature structure was obtained from (A23)-(A25) 
and is such that area C + area D = area E. 
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d d 

Normalized Forcing 



47 



Radius [km] 



Z= 0 




48 



V H0UR 12 




49 



U H0UR 12 




50 



f 

N H0UR 12 




H 

000 



L I 
-.079 - 



51 



T H0UR 12 




52 



V , H0UR 24 





53 



U H0UR 24 




54 



W H0UR 24 




55 





T H0UR 24 




56 



DELTA T 



H0UR 12 
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DELTR T 



H0UR 24 
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o 



0 

CJ 

1 



0 

1 



StiJi— ujq: cn 



0 

CD 

1 
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